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We study relaxation times, also called mixing times, of quantum many-body systems described 
by a Lindblad master equation. We in particular study the scaling of the spectral gap with the 
system length, the so-called dynamical exponent, identifying a number of transitions in the scaling. 
For systems with bulk dissipation we generically observe different scaling for small and for strong 
dissipation strength, with a critical transition strength going to zero in the thermodynamic limit. We 
also study a related phase transition in the largest decay mode. For systems with only boundary 
dissipation we show a generic bound that the gap can not be larger than ~ 1/L. In integrable 
systems with boundary dissipation one typically observes scaling ~ 1 /L 3 , while in chaotic ones one 
can have faster relaxation with the gap scaling as ~ 1/L and thus saturating the generic bound. 
We also observe transition from exponential to algebraic gap in systems with localized modes. 


I. INTRODUCTION 

With advancing quantum technologies [T| it is becom¬ 
ing increasingly important to understand interaction of 
quantum systems with external degrees of freedom. Evo¬ 
lution of a system coupled to environment can be de¬ 
scribed by a master equation. A particularly appeal¬ 
ing type of a master equation is a Lindblad equation - 
a rather general setting that can describe any Marko¬ 
vian evolution 2]. While Lindblad equations have been 
used extensively in the past to describe few-particle sys¬ 
tems in NMR or quantum optics, recently increasing ef¬ 
forts are devoted to understand many-body systems in 
the Lindblad setting. Motivation comes from a wide 
range of fields where Lindblad equations find their ap¬ 
plication, to name a few: as a computational resource 
in quantum information, to study e.g. transport proper¬ 
ties of strongly-correlated condensed-matter systems, or 
to study nonequilibrium statistical physics of many-body 
systems. 

Usually the object of most interest for open systems is 
a steady state, that is, a state to which any initial state 
converges after a long time. Besides the steady state, dy¬ 
namics is also of interest with one of the most important 
quantities being the relaxation time. In a finite system 
relaxation time is simply equal to the inverse gap of the 
Liouvillian propagator generating evolution. In statis¬ 
tical physics the scaling power of the gap is called the 
dynamical exponent - a critical exponent determining a 
universality class to which a model belongs. Dynamical 
exponents have been extensively studied in classical ex¬ 
clusion models, see e.g. Ref. [3! or Refs. 12 E] for some 
more recent results. 

In a quantum domain much less is known about re¬ 
laxation times of many-body systems. Depending on a 
situation one might want the gap to be large or small. For 
instance, if dissipation is engineered in order to prepare 
a specific steady state one might want relaxation to be as 
fast as possible. On the other hand, if dissipation is un¬ 
wanted, for instance, in a quantum memory device, relax¬ 
ation should be as slow as possible. The value of the gap g 
is important not just for the relaxation time itself 0 □, 


but can also carry information about the steady-state 
properties. Namely, if the gap is finite (so-called rapidly 
mixing systems) one can show that this implies a clus¬ 
tering of correlations in the steady state m , meaning 
that local observables are uncorrelated on a scale larger 
than ~ 1 /g. Rapid mixing also implies the stability of 
steady state to local perturbations fIT)f - lT2] . If the gap on 
the other hand closes in the thermodynamic limit this can 
lead to a nonequilibrium phase transition HMH] and can 
result in a non-exponential relaxation mm towards a 
steady state. Understanding how the gap scales with the 
system size is therefore of fundamental importance. 

There have been few scattered results in the litera¬ 
ture calculating or bounding the gap either analytically 
or numerically for specific Lindblad models. With our 
work we plan to extend these results, studying in more 
detail how the gap scales with the system size. One can 
distinguish grossly two different situations depending on 
the number of sites at which dissipation acts: in a lattice 
with L sites dissipation acts on ~ L sites (i.e., on most 
sites), a setting we will call bulk dissipation, or, it can 
act only on a fixed number of sites (number of sites does 
not grow with L), a setting we will call boundary dissipa¬ 
tion (eventhough sites at which it acts need not be at a 
boundary). Not surprisingly, as we shall see the gap can 
behave differently in the two cases. What is known so far 
is that for so-called Davies generators, also called ther¬ 
mal reservoirs, that are an example of bulk dissipation, 
one is in certain cases able to rigorously prove that the 
gap is independent of the system length EH 122 ],g~L°. 
One can also show exponential relaxation in weakly cou¬ 
pled systems [23] ■ Similarly, one can show that the gap 
can be constant for some other systems with bulk dissi¬ 
pation nzmsi, while in other systems it can also scale 
as ~ 1/L 2 H2HS1122- O n the other hand, for open sys¬ 
tems with boundary dissipation the observed gaps have 
so-far all been scaling as ~ 1/L 3 , or smaller, examples 
being the XY [201122 [25] or the XXX [27 ] [23] model. For 
scaling of gaps in the Redfield equation see Ref. [la¬ 
in the present work we are going to study the scaling 
of the gap with system size in a number of spin chain 
models with bulk as well as with boundary dissipation. 
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The aim is to get a better overview of possible gap scaling 
and changes in the scaling power as one varies parame¬ 
ters. Such changes can be associated with possible phase 
transitions in the steady state as well as in decay modes. 
In the main part of the paper we shall organize sections 
according to different models studied, explaining for each 
different techniques used (ranging from analytical to nu¬ 
merical) to infer the scaling. We shall also pay attention 
to the validity of a weak-dissipation perturbation theory 
that can be used to calculate the gap, demonstrating its 
failure in a number of cases. If one is not interested in all 
the details of gap calculation, or just wants an overview 
of different gap scalings found, there are Tables 0 and[IT| 
provided in the Summary section at the end of the paper. 

The paper is organized as follows. In Sec.II. we briefly 
explain the setting of Lindblad equations. In Sec III. we 
then study systems with bulk dissipation, while in Sec. 
IV. we study systems with boundary dissipation. Each 
of these two sections is split into subsections describing 
different systems. In Sec.IV. we also explain an argument 
that the gap in systems with boundary dissipation can 
not be larger than ~ 1/L. Finally, in Sec.V. we conclude 
as well as present a summary of the results found. 

II. LINDBLAD EQUATION 

The Lindblad equation is [291, 30; 

^=C(p):=i[p,H}+C dis (p), (1) 

where £dis(p) = %LjpL^ — pL^Lj — L^Ljp is a linear 
superoperator called a dissipator that can be expressed 
in terms of traceless Lindblad operators Lj. Denoting 
eigenvalues of the Liouvillian C by A j{£), and ordering 
them according to their real parts, with Ao(£) = 0 assum¬ 
ing to be nondegenerate, the gap is equal to a negative 
real part of the 2 nd largest eigenvalue, 

9 '■= — lRe[A 1 (,C)] (2) 

Eigenvectors of C corresponding to nonzero eigenvalues 
are called decay modes. 

The models that we are going to study will all be spin 
chains composed of L lattice sites, each carrying a spin- 
1 /2 particle. Coupling in H will always be local between 
nearest-neighbor sites only, and with open boundary con¬ 
ditions. Dissipation will also be local, i.e., each Lindblad 
operator Lj will act nontrivially only on a single site (dif¬ 
ferent Lj though can act on different sites). Two types 
of dissipation will be employed, both physically moti¬ 
vated. The first one is dephasing for which Lj ex crj, and 
which tries to destroy off-diagonal matrix elements (in 
the eigenbasis of cr z ). The second one is magnetization 
driving in which Lindblad operators proportional to rais¬ 
ing and lowering operators try to impose an imbalance 
in populations of spin-up and spin-down states. 

We mention that the many-body spin chain models 
studied here are within reach of present day cold-atom 


technology, e.g. Refs. GEQI32], with individual compo¬ 
nents like few qubit controlled dissipation [33] or Heisen¬ 
berg spin chains [331 E5] already demonstrated. 


III. BULK DISSIPATION 

A canonical model that we shall study in both bulk- 
and boundary-driven cases is the anisotropic Heisenberg 
model (XXZ model for short). For zero anisotropy the 
XXZ model goes into the XX model, which is especially 
simple and even allows for an exact asymptotic solution. 


A. XX with dephasing 

The Hamiltonian of the XX model is 

H = Y. a J a i+i+ a d a i+i- ( 3 ) 

3 = 1 

Dissipation £dis is given by dephasing of strength 7 act¬ 
ing independently on each site, which is described by a 
set of L Lindblad operators, 



We note that such Liouvillian C conserves the total 
magnetization Z = )>V cr z , that is, if | ip) is an eigen¬ 
state of Z with eigenvalue and similarly \ip) is an 
eigenstate of Z with an eigenvalue Z^, then C(\if>)(ip\) = 
^ 2 jk c 3 k\ 1 p 3 )( < Pk\ is a superposition of terms in which all 
tjjj and all ipk are again eigenstates of Z with eigenval¬ 
ues Z^ and Zy, respectively. More formally conservation 
means that UC(p)U i = £(UpW), with U being rotation 
around the 2 -axis. C therefore has a block structure and 
we shall label each block by a magnetization difference z 
and by a number of flipped spins r, 

z := Z^- Z v , r := (L- Z^)/2. (5) 

For an L-site chain the allowed values of 2 are from —2 L 
to +2 L in steps of 2, while that of r are r = 0,1,..., L. 
Of special interest will be sectors with 2 = 0 because 
they carry a steady state, i.e., an eigenstate of C with 
eigenvalue 0. A subspace with 2 = 0 and some value 
of r will be simply called an r-particle sector and has 

(operator) dimension (^) . Two most important ones 
are r = 1 (1-particle sector), being the smallest nontrivial 
one, and the largest one with r = L/ 2 (or r = (L ± l)/2 
for odd L) being called a half-filling sector because half 
of the spins are pointing up and half are pointing down. 

An important property of dissipative systems where H 
is quadratic in fermionic operators (via Jordan-Wigner 
transformation) and Lindblad operators are Hermitian 
(but not necessarily quadratic), and our XX model is an 
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example of such a system, is that equations for corre¬ 
lation functions split into a hierarchy of equations ac¬ 
cording to their order, i.e., the number of fermionic op¬ 
erators in the expectation value This enables 

one to exactly calculate few-point expectation values in 
the steady state, for instance in the presence of an addi¬ 
tional boundary driving [HEM, an incoherent bulk 
hopping mm, or special engineered dissipation m- 

Here we are not interested in the steady state but in¬ 
stead in the Liouvillian gap, in particular in its scaling 
with the system size L. Due to magnetization conserva¬ 
tion there are L + 1 steady states, one in each invariant 
subspace with z = 0 and a particular r. One can also 
easily see that such a steady state is a uniform mixture 
off projectors to all (^) diagonal basis states in the cor¬ 
responding r-particle subspace, for instance, for L = 3 in 
a 1-particle sector (r = 1, z = 0; subspace dimension 3 2 ) 
the steady state is ~ | 001)(0011 + | 010 )( 010 | + 1100 )( 100 |. 
For more details about such states, for instance their 
Schmidt spectrum, see Appendix [A] . Subspaces with 
z / 0 do not contain any zero eigenvalues (i.e., steady 
states) because they are orthogonal to the identity op¬ 
erator which is, due to trace preservation, always a left 
eigenvector of £ with eigenvalue zero and is therefore 
non-orthogonal to all steady states. Because we are in¬ 
terested in the gap we can limit our discussion to sub¬ 
spaces with z = 0 as they are the only ones that contain 
steady states. Subspaces with r = 0 and r = L are of 
dimension 1 and contain only the steady state and are 
of no interest to us. All other L — 1 subspaces contain 
a nontrivial gap. We shall be in particular interested in 
the smallest gap out of those, i.e., the global one, giving 
the slowest relaxation rate in the system. 

Because of a hierarchy of correlations we know that 
the eigenvalues, i.e. decay rates, of p-point correlation 
equations are equal to (some) eigenvalues of the Liouvil¬ 
lian Wl however, one in general does not know whether 
they also give the gap. Regardless of that they can be 
used as an upper bound on the global gap. For instance, 
for a model with an incoherent hopping the relaxation 
rate of 2-point correlations scales as ~ 1/L 2 [5B| mean¬ 
ing that the global gap can not be larger. Similarly, for 
a boundary-driven system numerical calculation of the 
global gap also resulted in ~ 1/L 2 scaling ITT. . 

By numerically diagonalizing our £ in each /'-particle 
sector we observe that the gap is in fact the same in all 
sectors. This in particular means that to calculate the 
global gap it is enough to consider the 1 -particle sec¬ 
tor which is of dimension (^) = L 2 (or its symmetric 

(L — l)-particle partner that has exactly the same eigen¬ 
values). This is due to a free nature of H 1 causing that 
the spectrum of the 1 -particle sector to be contained in 
higher-r sectors. One could be tempted to conclude that 
the 1 -particle sector will be rather trivial - this is cer¬ 
tainly true for a 1-particle sector of dimension L in the 
Hilbert space of states - however, here we are dealing 
with a 1-particle sector of dimension L 2 in the Hilbert 
space of operators. As we shall see this allows for a rich 


behavior, among other things for a discontinuous transi¬ 
tion in the scaling of the gap from a constant 1/L° for 
small dephasing strength 7 to ~ 1/L 2 for non-small 7 . 

Before going to the actual calculation of the gap let 
us pause for a moment and have a look at an alternative 
formulation of the eigenvalue problem for the whole £. 
The Liouvillian that acts on a (2 L ) 2 dimensional opera¬ 
tor space is a non-Hermitian linear operator that can be 
written as a “Hamiltonian” of a non-Hermitian spin-1/2 
ladder composed of L rungs. Each rung, spanned by 4 
pure states, takes care of one site of £ which is also 4 
dimensional (e.g., three Pauli matrices plus an identity). 
The resulting ladder is for the XX chain with dephasing 
very simple: it is composed of two XX-coupled chains 
along the ladder legs with an additional imaginary cou¬ 
pling along rungs, see Fig. [l] and Appendix of Ref. m 
for the mapping details. The steady state is an eigenvec- 

-Hxx 


—i(£ + 7 LI) = 


nxx 

FIG. 1. The Liouvillian superoperator for the XX chain with 
dephasing is equivalent to a non-Hermitian ladder Hamilto¬ 
nian, with XX interaction along upper and lower legs (thin 
lines) and an imaginary z — z interaction due to dephasing 
along rungs (double lines). 



tor corresponding to the eigenvalue zero, in other words 
a “ground state” of a non-Hermitian ladder. Because of 
dissipative coupling along rungs this ground state is al¬ 
ways, regardless of the value of 7 , a direct product of 
singlet states at each rung. The steady state is there¬ 
fore always dominated by dephasing 7 , forcing the steady 
state to be a singlet state. This is due to a very special 
structure of £ - no interaction in H (which plays an inert 
role) and a dephasing that kills all off-diagonal matrix el¬ 
ements. As we shall show, things are very different for 
the first decay mode determining the gap. There is a 
transition from a 7 -dominated phase to a different phase 
in which H becomes important. 

In the ladder formulation of £ we can also see why 
the 1 -particle superoperator sector is nontrivial: it corre¬ 
sponds (in an appropriate basis) to states with one par¬ 
ticle in the upper leg (bra) and one particle in the lower 
leg (ket). Therefore, a 1-particle superoperator problem 
is like a problem of 2 interacting particles on a ladder. It 
constitutes the simplest case of an interacting system. 


1. One-particle sector 

Let us now calculate the gap in the 1-particle sec¬ 
tor. Essentially the same eigenvalue problem, apart from 
different boundary conditions [55J, has been rigorously 
solved in Ref. |43j . Here we shall present a different and 
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FIG. 2. The gap g of the XX model with dephasing 7 = 
0.5. Shown is the exact gap (symbols), the full curve is the 
approximate formula j 8 j, while the dotted curve shows the 
asymptotic expression 27r 2 /( 7 L 2 ). The approximate result 
Q is accurate for L larger than the transition point L c := 
ivy/2/'y. 


approximate calculation of the gap, which is simpler but 
nevertheless accurate in the thermodynamic limit L —> 00 
in which we are especially interested. 

The basic idea is the following. We have seen that the 
steady state is always dominated by the dephasing and 
one can expect that in some range of dephasing strengths 
the first decay mode will also be of the same nature - that 
is dominated by 7 . If this is the case, then its eigenoper- 
ator will be close to diagonal because dephasing kills all 
off-diagonal elements. In the lowest approximation the 
unitary part £h couples diagonal elements \j)(j\, where 
| j) denotes a state of L spins with the j-tli spin being 
flipped to 11 ) while all others are in state 10 ), to two off- 
diagonal |j)(j + 1| and | j + 1)(j |. Such a tri-diagonal ap¬ 
proximation reduces the size of C in the 1-particle sector 
from L 2 to 3 L — 2. In addition, we can immediately see 
that on this three-diagonal subspace all (L — 1 ) states of 
the form \j)(j + 1 | + \j + l)(j| are eigenstates with eigen¬ 
value —4y (they are eigenstates of Ch with eigenvalue 0 
because of different signs when H acts from the left and 
from the right, while £dis(| 0 )(l|j) = — 2 ^ 10 )( 11 ^ and sim¬ 
ilarly for |1)(0| •). Reduction of C to (2L — 1) dimensional 
subspace spanned by | j)(j\ and | j)(j + 1 | — \j + 1 ) (j | has 
a block structure of form, 

/ 0 C T \ 

Tired = I 4 j) > (- 7 , k ~ 2 i y/^^Sjk 

( 6 ) 

Basis ordering is such that the first diagonal block cor¬ 
responds to L states \j)(j\, and the second to (L — 1) 
states | j)(j + 1| — \j + 1)0’|- Writing the eigenvector as 
(x, y) and the eigenvalues as A, and eliminating x, we get 
an eigenvalue equation CC T y = (A 2 + 47A)y, with C C T 
being a tridiagonal matrix of size (L — 1) x (L — 1) with 


matrix elements 

(CC T )j t j~ = 8(—2 dj : k + dj+i,k + (7) 

CC T is nothing but a standard matrix appearing is a 
solution of harmonic oscillators or a tight-binding model, 
and has eigenvalues —32 cos 2 ( nj/(2L )), j = 1,... L— 1, 
leading to eigenvalues A j of £ re d being A j = — 2y(l ± 
8 cos 2 ( 7 rj/( 2 L))/ 7 2 ). The gap is determined by 
the largest eigenvalue A j= l_i, and is 

9 = 2 7 ( 1 -y i -T sin 7dL ) ) s ^_ + ,... (8) 

The gap g gives the distance of the largest decay mode 
eigenvalue from the origin along a real axis, with the next 
largest decay mode being asymptotically at a distance 
4 g. For H with periodic boundary conditions the gap is 
4 times larger than the above g ([ 8 ]). We can see that for 
sufficiently small 7 the expression under the square root 
in Eq.([ 8 ]) can become negative. For large L this happens 
for 7 < 7 c, where the critical dephasing is 


Our approximate expression for the gap (§ is accurate 
only for 7 > 7 C , where the decay mode is indeed gov¬ 
erned by dephasing, see Fig. [2]. Note though that in the 
thermodynamic limit y c —> 0 and therefore Eq. (| 8 | be¬ 
comes exact for any 7 . One can see already in Fig.j2]that 
below 7 C (or L c , depending on which parameter is held 
fixed) the gap becomes independent of L and in the ther¬ 
modynamic limit equal to 4y (in this regime the largest 
decay mode eigenvalue is complex and the gap g is equal 
to its real part). Therefore, for the largest decay mode 
there are two phases: for 7 <C y c the decay mode is gov¬ 
erned by H with the gap being g = 4y, while for 7 > 7 C 
the decay mode is governed by dephasing and the gap 
is 3 x Tfgjs- One can use perturbation theory to show 
that the gap is indeed independent of system size for suf¬ 
ficiently small dephasing HU. The convergence radius of 
weak-coupling perturbation series for small 7 is 7 C © 
and algebraically shrinks to zero in the thermodynamic 
limit. Also, the thermodynamic limit L —> 00 and the 
weak coupling limit 7 —> 0 do no commute, which can be 
also seen in Fig. [3] showing the phase diagram. 

Finally, let us briefly comment also on the overall spec¬ 
trum of C. Numerically diagonalizing C on the 1-particle 
subspace the following picture emerges, see Fig. [4]. Most 
of the eigenvalues, ~ L 2 in number, are within the bulk 
laying in a complex plane. As one increases L eigenval¬ 
ues “evaporate” from the bulk and join a bunch of real 
eigenvalues to the right of the bulk. The number of these 
separated real eigenvalues is proportional to L. Decay 
modes corresponding to the separated real eigenvalues 
are dictated by dephasing, whereas decay inodes in the 
bulk are instead dictated by the Hamiltonian. Remember 
that in the absence of dephasing the spectrum of C would 
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FIG. 4. (Color online) Complex eigenvalues of £ for XX model 
with bulk dephasing, 7 = 1, A = 50, and 1-particle sector. 
For 7 > 7 c there are of order ~ A eigenvalues that are on the 
real axis and are separated from the bulk consisting of the 
remaining ~ A 2 complex eigenvalues. Real parts of the bulk 
eigenvalues are around —47 with the width being ~ 7 /A. 


FIG. 3. (Color online) Phase diagram of the largest decay 
mode, (a) Dependence of gap g /7 on dephasing and A. Full 
black line is the critical 7 C delimiting two phases. In (b) 
and (c) are shown expansion coefficients \ck.j | of the largest 
decay mode x, x = =1 c fc,ll^) 0 'l, showing an almost di¬ 

agonal dephasing-dominated phase with g ~ 1/A 2 in (b), and 
a delocalized phase with g ~ 1/A° in (c). 


be entirely on the imaginary axis. Such distinct nature 
of two decay mode types is in turn reflected in two differ¬ 
ent phases of the largest decay mode, Fig. [3]. Indeed, if 
one decreases 7 at fixed A, the real eigenvalues get “ab¬ 
sorbed” in the bulk, with the last one disappearing at 
ss 7 c, at which point a transition happens in the scaling 
of the gap with A. Approximate values of real eigenval¬ 
ues can be obtained from the approximation with £ re d 
([ 8 ]), resulting in a scaling of the j- th largest eigenvalue 
of £ as A j ~ ( j/L ) 2 . In the thermodynamic limit real 
eigenvalues therefore cluster around the origin with their 
density there having a square-root singularity. It has 
been observed nn that if the gap g closes in the ther¬ 
modynamic limit there is a possibility for an algebraic 
relaxation instead of an exponential one occurring when 
g is finite. Such an algebraic decay can be explained j2D] 
by clustering of eigenvalues around 0. In our XX chain 
with dephasing this clustering is of the same kind as in 
boundary-driven free models studied in Ref. [20] and as a 
consequence, in the thermodynamic limit the relaxation 
will have a power-law form. 


B. XXZ with dephasing 

In the XX model the bulk Hamiltonian describes non¬ 
interacting particles. Choosing the XXZ Hamiltonian 
instead, the additional coupling in the 2 -direction rep¬ 
resents the interaction between fermions in the Jordan- 
Wigner picture. In this section we shall therefore con¬ 
sider the XXZ spin chain with bulk dephasing. The 
Hamiltonian of the XXZ chain is 

H = Y, a >r +1 + a ]+i + - ( 10 ) 

3 =1 

with A being an anisotropy parameter. Dissipation is 
the same dephasing at each site as used already for the 
XX model ©■ Magnetization is again conserved and 
there is one steady state in each sector with magnetiza¬ 
tion difference 2 = 0 and r particles, see discussion for 
the XX model with dephasing. Using a Jordan-Wigner 
transformation the anisotropy part Acr|cr| +1 has a form 
~ Arijrij + -[ and therefore represents interaction between 
nearest-neighbor fermions. Analytical solution for eigen¬ 
values of C is not possible anymore and we will have to 
resort to various perturbation approaches and numerical 
calculation. 

Numerically calculating the gap one sees that, as op¬ 
posed to the XX model, the gap is this time different in 
different sectors. In particular, the global, i.e., the small¬ 
est gap, is from the half-filling sector with r = A/2 (or 
r = (A ± l)/2 for odd A) and not from the 1-particle 
sector, r = 1. We shall nevertheless first discuss the 1- 
particle sector, where in the thermodynamic limit things 
are simple. 
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FIG. 5. The gap g in the 1-particle sector of the XXZ chain 
with dephasing, 7 = 1 . Symbols show results of numerical 
diagonalization for anisotropies A = 0.5 and 2.0, while the 
full line is 27r 2 /L 2 , giving the asymptotic gap, which is equal 
to the one in the XX chain with dephasing Q. 


FIG. 6 . The largest non-zero eigenvalue ci of matrix R Jl2|| in 
the XXZ chain with dephasing, determining the gap Eq. ( ]13| ). 
Full line is 1/L 0 ' 8 , suggesting the asymptotic dependence for 

A > 1. 


1. One-particle sector 

Because the 1-particle sector’s size is L 2 one can nu¬ 
merically calculate the gap for systems of reasonable size. 
Results are in Fig. [5]. As we can see for large system sizes 
the gap is equal to 


2tt 2 

7^2’ 


( 11 ) 


and is therefore equal to the one for the XX chain with 
dephasing. In the 1-particle sector the gap is in the ther¬ 
modynamic limit independent of interaction (anisotropy) 
A. This can be qualitatively understood if one looks at 
the non-Hermitian ladder formulation of £ (for the XXZ 
chain the only difference compared to the XX chain in 
Fig. [I] is that the Hamiltonian along the two legs is the 
XXZ Hamiltonian, see Ref. EH)- Because there are only 
two particles present in the ladder interaction is infre¬ 
quent, it is in fact only a ~ 1/L boundary effect (44], and 
is asymptoticall y irr elevant (for finite L correction to the 
asymptotic gap (11 1 due to finite A scales as ~ 1 /L 4 ). 


2. Global gap 

As mentioned, for the XXZ chain with dephasing the 
smallest gap is from the so-called half-filling sector with 
r = L/2 particles, which is also the largest subspace of 
C. Here we shall study the gap in this sector. We are 
going to use perturbation theory for small 7 as well as 
for large 7 , while in-between numerical calculations will 
be used to infer a general form of g. We are also going 
to show that the perturbation theory for large A fails. 

Let us start with small 7 . For the unperturbed part 
of the Liouvillian we take the whole unitary part, £q := 


i [p,H\, generated by the XXZ Hamiltonian (10), while 
perturbation is the dephasing £1 := £dis- Because the 
unperturbed £0 has a degenerate steady state subspace 
corresponding to eigenvalue 0 , we have to use 1 st or¬ 
der degenerate perturbation theory. The steady-state 
subspace is spanned by all projectors Xk = iV’fcXV’fcl 
to eigenstates | ipk) of H (we assume a non-degenerate 
spectrum of H). Let us denote a projection of £1 to 
the steady-state subspace of £ 0 by 7 R := VC\V, with 
Rj,k = tr(xj£i(xfc)) being of size (^). Using the form of 
dephasing Lindblad operators Q we get 


11j . k — LS 


1 j,k 


■£l<Vb 

p =1 


* z P m\■ 


( 12 ) 


The largest eigenvalue of R is equal to 0 and the next- 
largest one, denoted by —c\ (similarly as for £, all eigen¬ 
values of R have non-positive real parts), determines the 
gap of £ for small 7 , 

g = 701 + 0 ( 7 2 )- (13) 


We observe that, while on one hand dephasing can be 
derived as being due to a classical fluctuating magnetic 
field in the 2 -direction [55], we also see that the gap of 
£ for small dephasing is determined by eigenstate fluc¬ 
tuations of magnetization in the 2 -direction, as reflected 
in the form of R (12). In Fig. [d] we show results of exact 
numerical diagonalization of R for different anisotropies 
and system sizes, all in the largest half-filling sector with 
r = L/2. We can see that for A < 1 the eigenvalue C\ be¬ 
comes asymptotically independent of £, for instance, for 
A = 0 it converges towards 4 as already calculated in sec¬ 
tion about the XX model. For A > 1 it on the other hand 
decays; numerical data for available L < 18 is consistent 
with a ~ 1/L 0 8 decay. For small (A — 1) this asymp¬ 
totic decay starts for sizes larger than L* oc 1/(A — 1). 
It would be interesting if one would be able to directly 
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FIG. 7. (Color online) The largest non-zero eigenvalue in the 
even/odd sector for the XXZ chain with dephasing, all for L = 
8 . Insets show dependence of the critical -y c that determines 
validity of perturbation series. Left: gapless regime of A = 0.5 
for which 7 C depends on L algebraically (dashed curve in the 
inset is ~ 1/L 1 ' 2 ). Right: gapped regime A = 2.0 for which 
7 c x exp (— kL ) (inset). 


calculate the largest nontrivial eigenvalue of R via the 
Bethe ansatz eigenstates 'i/’fc- 

So far we have calculated the global gap for small de¬ 
phasing 7 without saying anything about the convergence 
radius, i.e., the validity of such an approximation. To find 
out critical y c upto which perturbative gap (13) can be 
expected to hold, we are going to numerically calculate 
the gap and compare it to perturbative prediction. In ad¬ 
dition to symmetries of £ already pointed out there is also 
a spatial reflection symmetry with respect to exchanging 
sites j —> L — j, resulting in even and odd decay modes 
and eigenvalues of £. The largest nonzero eigenvalue 
in even and odd sector (both for 2 = 0 and half-filling 
r = L/2) is shown in Fig. [7] for system size L = 8 and 
two values of A. The steady state - a uniform combina¬ 
tion of all pure-state projectors - is always from the even 
sector. We can see that the largest non-zero eigenvalues 
cross at a critical y c . The decay mode that determines 
the gap in the half-filling sector (red pluses in Fig. [7]) is 
from the even sector for 7 < y c and from the odd one for 
7 > 7 C . For 7 < 7 C we see that the gap (red pluses) is 
proportional to 7 (horizontal line in Fig. [7]) and therefore 
perturbation series (13) holds. The convergence radius 7 C 
shrinks algebraically with L in the gapless phase A < 1, 
at A = 0.5 it decays as y c x 4.1/L 1 - 2 , although we can 
not exclude asymptotic ~ 1/L scaling, while it is expo¬ 
nentially small in the gapped phase A > 1. For the XXZ 
chain with bulk dephasing the convergence radius of per¬ 
turbation series in the dephasing therefore shrinks to zero 
in the thermodynamic limit. We also note that the two 
eigenvalues shown in Fig. [7] are both real, except the odd 
one for 7 < y c and A = 0.5 which forms a complex pair. 


For non-small 7 the XXZ chain with bulk dephasing is 
diffusive ) 46 | and therefore one in general expects the gap 
to scale as ~ 1 /L 2 in the thermodynamic limit. Namely, 
for diffusive systems evolution of macroscopic observables 
does not depend independently on time t and spatial co¬ 
ordinate x but instead only on one scaling variable x 2 /t 
and therefore time should scale with a size squared. We 
numerically calculated m the gap for fixed non-small 7 



FIG. 8 . (Color online) Gap for the XXZ chain with dephasing 
in the half-filling sector, 7 = 1. Full lines denote asymptotic 
~ 1/L 2 scaling for all A. 


and indeed confirmed the g ~ 1 /L 2 scaling, see Fig. [8] as 
well as Ref. [ 19 ]. Exact dependence of g on 7 and A is 
more complicated; a crude approximation that seems to 
work for large L (at least for 7 = 1 ) is g ~ 4 7 2+a 2 “ 

the straight lines plotted in Fig. [H] is in fact this expres¬ 
sion. 

When some parameters in the system are large one 
can again use perturbation theory. One such possibility 
is the case of large dephasing 7, where the unperturbed 
Liouvillian £0 contains Hq = Yj ^■ a j a j +1 an d dephas¬ 
ing dissipation £4^ while the perturbation C\ consists of 
hopping given by Hi = Yj ( 7 j cr j+i + rftf+i- This has 
been done in Ref. [ 19 ], obtaining that the gap for large 7 
scales as ~ 1/(7L 2 ). 

We shall redo calculations of Ref. mi deriving also 
subleading terms, in order to be able to comment on 
the validity of different limits. First, £0 has a degener¬ 
ate steady-state subspace and one has to use degener¬ 
ate perturbation theory. First-order perturbation on the 
steady-state subspace of dimension 2 L is always zero be¬ 
cause all steady states are diagonal, while the hopping 
£1 transports one spin in the bra or ket, and so at least 
two applications of £1 are needed to again get a diago¬ 
nal state and with it a nonzero matrix elements of C\. 
Second-order perturbation of the steady state manifold 
is determined by, 


£ e ff = —'P£ 1 £q 1 £ 1 'P, (14) 

where V is a projection operator to the steady state man¬ 
ifold. The steady state manifold of £ 0 consists of all 
diagonal density matrices, \ipj)(ipj\, where \ipj), j = 
1 ,... 2 l are 2 L (basis) states in the standard er z - 
eigenbasis. One can write matrix elements of £ e ff as 
[H e s]j,k ■= where matrix 

Heff can in turn be written in terms of Pauli spin vari- 
























ables, obtaining 
Hxxx , 


H eS = 


+ 


7 

A 2 /7 


A 2 /7 
A 2 + 47" 


[(Tj • (72 + <JL -1 ' — 2 • 1] + 


L -3 


2 (A 2 + 7 2 ) ^ 


- 1)(1 - 07+1 • 07+2)5(15) 


where Hxxx := X)y=i(l ~ <?j ' 07 + 1 )- Such H eS is Her- 
mitian, has L + 1 zero eigenvalues (one for each invari¬ 
ant magnetization sector), while all other eigenvalues are 
positive. The smallest non-zero eigenvalue is equal to 
the gap g of £ in the limit when C\ -C Cq. In Ref. nu 
the above result (15) has been derived in the leading or¬ 
der in I/7, i.e., H e g = Hxxx/ 7> from which due to a 
quadratic low-energy dispersion of Hxxx one gets the 
asymptotic gap g x -Jp- This perturbative result is also 
valid in the thermodynamic limit m because the con¬ 
vergence radius does not shrink to zero. With the exact 
expression for H e g we can also explore the case of large 
A, for which one would be tempted to think that pertur¬ 
bative series will again work because C\ is again small 
compared to £0 that contains a large parameter A. Such 
reasoning though is in fact wrong and one can not use 
perturbation in C\ if only A is large. The reason for the 
failure is rather instructive and we are going to explain 
why it happens. As opposed to the limit 7 — > 00, for 
A — > 00 the two terms that appear in addition to Hxxx 
in Eq. (15) can not be neglected because the two pref¬ 
actors in front of them scale for large A as I/7 and are 
therefore of the same order as the Hxxx term. There 
is also an additional subtlety: it would be tempting to 
retain only the leading order expansion of the two pref¬ 
actors, the already mentioned I/7 and 1/(27), however, 
in that case the ground state of H e g is highly degenerate 
and the gap would therefore remain zero (one of the rea¬ 
sons for the degeneracy is that the boundary term in the 
first line of Eq. ( fl5] ) with a I/7 prefactor exactly cancels 
the boundary terms from Hxxx, leaving the first and the 
last spin uncoupled). One has to retain at least the 1st 
order expansion of the prefactors, resulting in terms that 
are proportional to This means that the gap of H e g , 
and with it also the Liouvillian gap g , will scale as g oc 
for large A (this conclusion remains true despite the fail¬ 
ure of 2nd order perturbation expansion). Calculating 
the gap of the full H e g (15) and comparing it with the 
exact gap of £, one gets data in Fig. ■ As one can see 
the error does not decrease to zero as one increases A 
(this residual error as well increases with L). The rea¬ 
son is that for large A and fixed dissipation 7 not all 
non-zero eigenvalues of £0 are large - some are of order 
7 - and therefore a pseudoinverse £q 1 in perturbation 
series (14) is not necessarily small. This occurs because 
the spectrum of H 0 is highly degenerate with eigenstates 
being product states in the eigenbasis of a z , which is also 
the eigenbasis of £di s - If | ip) and \y) are two such degen¬ 
erate eigenstates we will have [\il>)(ip\, Hq\ = 0 while for 
dephasing ||£ dis (|'0)((/?|)|| ~ 7. As a consequence, con¬ 
secutive orders in perturbation series do not necessarily 



FIG. 9. (Color online) Relative error of the gap calculated 
from perturbative H e g (15 I, 7 = 1. For large A the error does 
not go to zero, signifying that the second-order perturbation 
theory in large A fails. 


uniformly decrease. This can be compared to perturba¬ 
tion series for fixed A and large 7, for which all non-zero 
eigenvalues of £0 are large and of order 7 and therefore 
£q 1 is small. One can in fact draw a general rule: if 
a large perturbation parameter is in global dissipation 
(strong coupling) a perturbative series will be well be¬ 
haved, whereas if a large parameter is only in Hq one 
must be careful if Hq has degeneracies. Similar compli¬ 
cations can occur if (strong) dissipation acts only on a 
few sites. 

Let us summarize our findings for the XXZ chain with 
bulk dephasing and the half-filling sector: for small de¬ 
phasing 7 < 7c the gap is ~ 1 /L° for A < 1, while 
it scales as ~ 1 /L 08 for A > 1 . Critical dephasing 7 C 
decays algebraically with L for A < 1 while it is expo¬ 
nentially small for A > 1 , see also Table [I]. Perturbation 
series in small dephasing therefore fails in the thermody¬ 
namic limit. For non-small dephasing the gap is ~ 1 /L 2 
regardless of A, as one would expect for a diffusive sys¬ 
tem. Perturbation series for large 7 works, while it fails 
if A is the only large parameter. 


C. Constant gap 

We have seen in the XX and XXZ models that the gap 
can be constant for a sufficiently weak bulk dissipation. 
Problem with these two cases is that the critical dissipa¬ 
tion goes to zero in the thermodynamic limit. A natural 
question is: can one have a constant gap also for non¬ 
small dissipation? The answer is yes and we are going to 
give a simple example. Known are examples with only 
dissipation and no Hamiltonian PH 05 ]. An example 
we are going to present has a nonzero Hamiltonian and 
nonzero dissipation. It is a XX chain with an incoherent 
“hopping” given by Lindblad operator 

L 0 = a t (T j+ 1 > 


( 16 ) 














9 



10 100 200 
L 


FIG. 10. (Color online) The gap g in the 1—particle s ecto r 
of the XX chain with an incoherent one-way hopping (161. 
The kink at L = 10 is because the eigenvalue responsible for 
g goes from being complex to real. 


at each site j, = (cr* ± icrJ)/2. The Liouvillian C 
again conserves magnetization difference 2 and particle 
number r. Numerically diagonalizing C in a r = 1 parti¬ 
cle sector (z = 0) one gets gaps shown in Fig. 10 . We can 
see that asymptotically the gap is independent of L. Be¬ 
cause in the 1-particle sector interaction asymptotically 
does not matter, the same asymptotic gap would be ob¬ 
tained also in the XXZ chain with an incoherent hopping. 
Beware that if we would take the XX chain with an inco¬ 
herent hopping in both directions (i.e., adding additional 
Lindblad operators 1 ) the gap in the 1-particle sec¬ 

tor would asymptotically decay as g x 4tt 2 /L 2 , see also 
Ref. j55j . 

Another way of having a constant gap is to use a so- 
called thermal dissipators (also known as Davies genera¬ 
tors [19]) that can be derived in the limit of weak coupling 
to thermal reservoirs. The stationary state of such a mas¬ 
ter equation is thermal, though with Lindblad operators 
that are not local. For such models one can in certain 
cases rigorously prove that the gap is constant |22j . 


IV. BOUNDARY DISSIPATION 

In this section we are going to study the gap in open 
systems with boundary dissipation, that is, with Ldis 
acting in the thermodynamic limit nontrivially only on 
a finite number of sites. In all our cases H will al¬ 
ways be a spin—1/2 chain and dissipation will act on 
the leftmost and rightmost lattice sites. As mentioned, 
for boundary-driven open systems the gaps observed in 
the literature are all ~ 1/L 3 or smaller. Examples 
are ~ 1/L 3 in the XY chain with boundary magneti¬ 
zation driving p~8l l20l [26] (or ~ 1/L 5 at nonequilibrium 
phase transition points [25]; or even ~ 1 /L 7 on the so- 
called resonances [ 50 ]); ~ 1/L 3 is the scaling also for the 
magnetization-driven XXX model [23], or for the XXZ 


model with an incoherent hopping as a driving ]2§) . 

Considering these results one can ask whether the gap 
in a boundary-driven open system is perhaps always ~ 
1/L 3 or smaller? As we will see the answer is no. But 
let us first present an argument that the gap can not be 
larger than ~ 1/L. 


A. Gap upper bound 


Let us have a one-dimensional system described by a 
local Liouvillian, C = X)y=i £j> which is a sum of local 
term Cj , each of which acts nontrivially only on a fixed 
number of consecutive sites around site j (Hilbert space 
dimension of a single site is finite, implying that Cj are 
bounded, and the steady state is assumed to be unique). 
In addition, let all but a fixed number of Cj be purely 
unitary, i.e., dissipation is present only in a fixed num¬ 
ber of Cj. In the thermodynamic limit the number of 
consecutive sites without any dissipation therefore grows 
linearly with the system length L. The question we want 
to ask is: what is the fastest possible relaxation time, i.e., 
the largest possible gap, in such an open system? 

We will show that the gap can not be larger than 
~ 1/L, or, in other words, relaxation can not happen in 
a time that grows with L slower than linearly. The argu¬ 
ment is actually very simple. Because there are sites that 
are of distance ~ L away from the nearest site with dissi¬ 
pation, a “disturbance” at that site can not dissipate in a 
time smaller than oc L, i.e., the time a disturbance needs 
to get to a site with dissipation. We can also use a trans¬ 
port argument: for unitary evolution local conservation 
of energy holds and therefore, because the local energy 
current is a bounded operator, it will take at least a time 
oc L for the energy of the initial state to be dissipated, if 
we choose an initial state having total energy propor¬ 
tional to L. One could also rigorously formulate the 
above argument by e.g. using the Lieb-Robinson bound 
for open systems [SHU EE]. The Lieb-Robinson bound 
essentially formalizes a statement that there is a finite 
propagation speed in bounded locally-coupled systems. 
One consequence is that connected correlations of local 
operators get exponentially suppressed outside of a light 
cone, or, that the Heisenberg picture A(t) of the initial 
local operator can be approximated by a part of C with 
support inside the light cone (outside the light cone one 
has A(t) w 1). Taking a local A(0) J_ 1 with support on 
sites that are a distance ~ L away from dissipation, one 
immediately sees that, because one has A(t -A oo) = 0, 
relaxation time can not be smaller than ~ L. 

In next two subsections we are going to study by now 
familiar XX and XXZ models, showing that the Lieb- 
Robinson bound g ~ 1/L is never reached, and then in 
the last subsection show some examples for which one 
does get g ~ 1/L. 
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B. XX with boundary dephasing 

We take the XX chain (|3| with dephasing on the first 
and the last sites, that is, with only two Lindblad opera¬ 
tors L\ and Ll from Eq. Similarly as in the XX chain 
with dephasing on all sites, the gap is again the same in 
all r-particle sectors and we can limit our discussion to 
the 1 -particle sector. 


4 
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FIG. 11. Left: Eigenvalues of C in the 1-particle sector for 
the XX chain with boundary dephasing of strength 7 = 0.01. 
Three crosses are eigenvalues that have the largest real part 
and give g ( |18| ). Right: Numerical g( 7 ) (symbols) for L = 40 
and prediction by Eq.( 191 (full curve). 


We are first going to use perturbation theory for small 
7 . An example of a spectrum of C in the 1-particle sec¬ 
tor is in Fig. m- One observes that for small 7 the 
largest eigenvalue, determining the gap, is complex and 
has the largest absolute value of the imaginary part (two 
crosses in left Fig. 11). This pair of complex eigenval¬ 


ues determines the gap all the way upto the maximum of 
5 ( 7 ) (which for L = 40 shown in right frame of Fig. |TT] 
happens around 7 « 1). For larger 7 the gap is given 
by a real eigenvalue (the cross on a real line in Fig. © 
which though has in the leading order the same real part 
as the complex pair. This means that one can use sim¬ 
ple nondegenerate perturbation theory to get the gap for 
small 7 . The eigenvalue with the largest (or smallest) 
imaginary part of Ch is po = |' 0 i)(' 0 l| j where we denote 
eigenvalues of H by Ek = 4cos q k ,h = 1 where 

9k ■= 2 ( 77 , and eigenstates by \ijj k ) = J2j c jk\j), c r = 
yj2/{L + 1 ) sin q r , where | j) denotes a state with all spins 
down apart from the j-th spin. The eigenvalue correction 
will then be given by n := tr(po£di s (/ 0 o))) where dephas¬ 
ing £dis is a sum of nontrivial parts acting on the first 
and the last sites. £dis gives a non-zero value (equal 
to — 27 ) only when acting on a non-diagonal operator 
on the first/last site, and therefore £dis(Po) is a sum of 
terms that are either |j)(l| or \j)(L\ (or their Hermitian 
adjoint), i.e., oc — 27 ^-^ c\c*\l){j\ + • • •, at the end 
resulting in 

87 


(L+iyi 


- (2 - cos q 2 - cos g 2 i 2 ) sin 2 q L + 


L 

+ Sfe.L + Sfc, 1 + Sqfc + SL,k 
fc=1 


(17) 


where Sk, r := sin 2 q k sin 2 q r ^. For large L the term in the 
first line of Eq. is subdominant, while the second 


line gives a large-L expression for k and therefore also 
the gap, resulting in 


9 = 1 


167 T 2 


(L + iy 


0 ( 7 2 ). 


(18) 


It is instructive to qualitatively understand how the ~ 
1 /L 3 scaling comes about: we see that s k , r is just a prod¬ 
uct of eigenstate expansion coefficients Cj and that the 
sum scales as J2k l c i| 2 l c fci| 2 ~ |ci| 2 , which in turn is the 
smallest overlap probability at the first site, giving one 
1 /L due to normalization and an additional qf ~ 1 /L 2 
due to the longest-wavelength eigenstate which is the 
slowest to relax. Such a scenario has already been ob¬ 
served j2D] in the XX chain with boundary injection of 
particles, i.e., a model with Lindblad operators ~ cr ± 
at the boundaries instead of our ~ cr z , making it ex¬ 
actly solvable. One can in fact argue that such ~ 1 /L 3 
scaling is generic for boundary driven models with inte- 
grable Hamiltonian Hq having a plane-wave like longest- 
wavelength eigenstates. 

For large 7 one on the other hand expects the gap to 
decay with dephasing strength as ~ 1/7 because of an 
effective decoupling of the system and as a consequence 
increasing relaxation time. One can extend snrall -7 ex¬ 


pression (18) to also have the correct large -7 behavior by 
writing, 


9 = 


1 


16t r 2 


1 + 7 2 L 3 


(19) 


One can see in right Fig. EH that this expression indeed 
fits well numerical results (small discrepancy seen is due 
to a subleading ~ 1/L 4 correction to Eq. (19)). Observe 
that for the boundary driven XX system, as opposed to 
the XX with bulk dephasing, there is no transition in the 
scaling of the gap as one varies 7 . 


C. XXZ with boundary dephasing 

For the XXZ model with boundary dephasing we take a 
standard XXZ Hamiltonian (10) and the same two Lind¬ 
blad operators L\ and Ll Q. We begin by discussing 
weak dephasing. 


For weak dephasing one observes that the gap comes 
from real eigenvalue of C and one therefore has to use de¬ 
generate perturbation theory in an appropriate r-particle 
subspace (again, always limiting to z = 0). First order 
already gives a non-zero contribution and so the proce¬ 
dure is very similar to the one already used for the XXZ 
chain with bulk dephasing, leading to Eq. (12). Repeat¬ 


ing the same steps now for £dis that acts only on the 
boundary two sites, one gets 

Rj,k = \(^Ml\^k )\ 2 + \(yjWl\i’k )\ 2 - 2 5 jtk: , ( 20 ) 

where \ip k ) are eigenstates of the XXZ chain in an appro¬ 
priate r-particle sector. Eigenvalues of R are for small 7 
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FIG. 12. (Color online) 1-particle sector of the XXZ chain 
with boundary dephasing, (a) Numerically calculated gap g 
for small 7 = 0.01 and the theory (full curves, Eq. HD), (b) 
g( 7 ) always scales as ~ 1/E 3 ; symbols are numerical values 
for two different sizes and A, full curves are Eq. \22\. 


equal to the largest real eigenvalues (divided by 7 ) of C. 
Doing the calculation in the 1-particle sector one gets 


9 = 


16t r 2 


7 


(1 + E) 3 (1 + A)- 


+ 0 { 7 2 ). 


( 21 ) 


For small 7 and in the 1-particle sector the gap therefore 
scales as ~ 1/E 3 irrespective of the value of A. For a 
non-small 7 expression that fits numerics well and has a 
correct weak-depliasing limit ( 21 ) is 


167T 2 7 

(E + l ) 3 y 2 + (1 + A ) 2 ' 


( 22 ) 


We see in Fig. [T 2 Jd that the perturbative result in 
Eq. ©, i.e., horizontal dependence for small 7 in the 
figure, holds up-to an E-independent 7 . 

We now move to the half-filling sector. In Fig. [13| we 
show the results. For large E the gap scales as g ~ 1/E 3 
in the gapless phase A < 1, while it is exponentially 
small, g ~ exp (—aL), in the gapped phase of A > 1 [52] , 
Comparing gaps in the 1-particle and the half-filling sec¬ 
tor one also observes (data not shown) that for A < 1 
the gap in the 1 -particle sector is the smallest of the two, 
while for A > 1 the smallest gap is from the half-filling 
sector. For arbitrary 7 the global gap therefore scales as 
~ l/L 3 for A < 1 and ~ exp (—aL) for A > 1 , see also 
Table HU. 

An exponentially small gap can be most easily under¬ 
stood for small 7 via matrix R. Without loss of gen¬ 
erality we limit to even L and focus on the half-filling 
sector. For matrix elements of R eigenstates of H and 
matrix elements of Lindblad operators is what matters. 
In the spectrum of H there are two almost degenerate 
eigenstates 2 which are the most important for the 
gap, i.e., for the largest nontrivial eigenvalue of R. For 
large A they are a symmetric and antisymmetric com¬ 
bination of a domain wall in the middle of the chain, 
\ 1 jj 1 ) ~ | R) + | L) and (^ 2 ) ~ |-R) — | L), where we denoted 
| L) = |11...100...0) and | R) = |00... Oil... 1) (| L) is a 
state with the left-half of spins down, and |i?) a state with 
the right-most half of spins down). For finite A domain 
wall states |.R, L) differ from a perfect wall only within a 
localization length ~ 1 /lnA of the middle spin, see e.g. 
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FIG. 13. (Color online) Half-filling sector of the XXZ chain 
with boundary dephasing, (a) The largest non-zero eigen¬ 
value of R, Eq. (20), determining the gap for small 7 (131. 
Asymptotic behavior is denoted by dotted line (= 139/E 3 ) 
and full curve (~ exp (—1.32E)). (b) The gap for non-small 
7 = 1 . Asymptotically the scaling is the same as for small 7 , 
namely, the dotted line is 80/E 3 (A = 0.5) and the full curve 
is ~ exp (—1.32E) (A = 2.0). 


the Appendix in Ref. [53]. For A > 1 the two states 2 
are therefore “localized” around the site n/ 2, their energy 
is almost degenerate, they are the highest energy states 
(in the half-filling sector), and they are gapped by ss 2A 
from the rest of the spectrum. They have an important 
property that <jf\ijji) « |t/i 2 ) and ~ i.e., a 

subspace spanned by 7 / 7,2 is invariant under a\. As a 
consequence, overlaps (^*^ 1,2 Iff^ 1 , 2 ) are exponentially 
small in E and in the matrix R a 2 x 2 block correspond¬ 
ing to these two eigenmodes is decoupled from the rest. 
This 2x2 block is of the form —21 + (2 — e)cr x resulting 
in two eigenvalues —4 + e and — e, with the correspond¬ 
ing Liouvillian decay eigenmodes IV’i) (^i I + for 

eigenvalue —e, and |t/>i)(t/>i| — \'ip 2 )(ip 2 \ f° r a symmetric 
partner at —4 + e (remember that the steady state is a 
uniform mixture of all projectors, J2k \' l Pk)(' t l’k\)- The gap 
is thus equal to e which is in turn determined by the lo¬ 
calization length of two localized eigenmodes. The rate a 
in g ~ exp (—aL) is therefore proportional to the inverse 
localization length, resulting in a oc In A. Exponentially 
slow relaxation found for A > 1 could be interesting for 
instance for quantum memory - we see that the domain 
wall can support a one-qubit quantum memory formed 
out of 7 / 7,2 which is exponentially resilient to dephasing. 
Unfortunately, as we shall see in the next magnetization- 
driven model, this resilience is lost for boundary dissi¬ 
pation in a transverse direction. We note that one can 
also get exponentially small gap because of localization 
due to disorder, an example being a boundary driven XY 
chain with a disordered magnetic field [25] , 


D. Magnetization-driven XXZ 

So-far we have used dephasing as our canonical exam¬ 
ple of dissipation because it eased up theoretical anal¬ 
ysis. Steady states are in those cases rather simple, 
namely uniform mixtures of all projectors in a respec¬ 
tive symmetry subspace (Appendix [A|. Much more in¬ 
teresting steady states arise if one has boundary driving 
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FIG . 14. The gap g for the XXZ model and boundary driving 
(231 with /u = 0.1 and r = 1. Straight dashed line suggests 
asymptotic scaling g oc 1 /L :i . 


that breaks symmetries, for instance imposing a left-right 
asymmetry by having different driving at chain ends. In 
such a case the steady state will be a genuinely nonequi- 
librium state with nonzero currents flowing through the 
system. 

In this subsection we are going to study one such ex¬ 
ample, on one hand to contrast the obtained scaling with 
that of a boundary dephasing case and on the other hand 
to calculate the gap for a physically much relevant open 
system. We shall take the same XXZ model as in previ¬ 
ous sections |To| ), without dephasing, but instead with a 
boundary magnetization driving described by 2 Lindblad 
operators at each end, 

Li = y/r(l + n) af, L 2 = \/r(l - fj.) erf, 

L 3 = y/T(l-»)*+, L A =y/r(l + fi )aZ. (23) 


Driving parameter /r parametrizes the asymmetry be¬ 
tween injection and absorption of particles at the bound¬ 
ary, trying to impose expectation value of cr z equal to ±/x 
at chain ends, while T is the coupling strength. Nonzero /r 
therefore causes a nonzero magnetization gradient along 
the chain. We are going to use g, = 0.1, however, the 
scalings observed (e.g., Table 0 are independent of the 
precise value of // as long as g is not too close to g = 1 
for which one gets a blockade of transport due a step-like 
magnetization profile resulting also in an exponentially 
small gap leading to exponentially slow relaxation [54] . 
The Liouvillian of such magnetization-driven XXZ model 
still conserves z. but not anymore r ([5]). Note that T plays 
the role of dissipation strength, similarly to 7 in the case 
of dephasing dissipation - small T means weak dissipa¬ 
tion, e.g., weak external coupling. Dissipation strength in 
the paper is therefore determined either by 7 for dephas¬ 
ing or by T for “magnetization” driving in Eq. ( |23| ). 

Numerically calculated gaps are shown in Fig. [14] . The 
gap g looks to have a nice ~ 1/L 3 scaling regardless of 
the anisotropy A. We can compare the gap obtained 
here to the one obtained in previous subsection for the 


XXZ model driven with boundary dephasing (Fig. 131: 
in the gapless phase of A < 1 the scaling is in both 
cases ~ 1 /L 3 , while in the gapped phase of A > 1 the 
gap is here ~ 1/L 3 while it was exponentially small for 
the boundary dephasing case. We see that the scaling of 
the gap can change already by changing boundary terms 
only. 

Analyzing the gap for small coupling T, i.e., small dis¬ 
sipation, again requires application of a degenerate per¬ 
turbation theory. Using Lindblad operators (231 the per¬ 
turbation matrix R is in this case 


R jk = 2(1 + ^)|(^'|crj h |V’fe )| 2 + 2(1 - n)\{ipj\a 1 | ip k )\ 2 + 

+ 2(1 - n)\tyjWi IV’fe)! 2 + 2(1 + M)l(^>Zl^>| 2 - 

- 2 fi 6 jk ('ipj\(rl - a z L \^3) ~ 4 < 5 jfc- ( 24 ) 

The spectrum of such R for the XXZ chain is in fact 
independent of g. For g = 0 the above R can be further 
simplified, taking into account also that H is real, results 
in 


Rjk = 2\(ipj\al\ipk)\ 2 + ^\{^jWl\^k)\ 2 - 4Ay fc . (25) 

Such R is real and symmetric with the eigenvector corre¬ 
sponding to eigenvalue 0 being a uniform superposition 
of all basis states (we remind that perturbative R has 
always one eigenvalue equal to 0 with the corresponding 
eigenvector though being in general more complicated). 
The gap (i.e., C\) for small coupling T of magnetization 
driving therefore depends on fluctuations of <r x at the 
boundary two sites; this can be contrasted with the de¬ 
phasing case (20) where fluctuations of cr z matter. In 



FIG. 15. (Color online) The XXZ chain with magnetization 
driving and small dissipation T. (a) The largest non-zero 
eigenvalue ci of R (Eq. ( |24[ )) determining the gap for small F. 
Full line suggests asymptotic scaling ~ 1/L 3 . (b) Dependence 
of exact gap on T and perturbative ci from frame (a), all for 
A = 1.5 (for A = 0.5 behavior is similar). Deviations from c 1 
begin at an L-independent T. 

Fig. [15] we show results for c\ and verification of its va¬ 
lidity. We can see that for small T the gap also scales as 
~ 1 /L 3 irrespective of A and that approximation g ss Tci 
holds upto an L-independent T. What is interesting is 
that, compared to the XXZ model with boundary de¬ 
phasing we changed only boundary driving, that is, we 
modified only action of £ on 2 out of L sites, and never¬ 
theless the scaling of gap for A > 1 changes from expo¬ 
nential to algebraic. On the level of perturbative matrix 
R this can be understood as being due to the breaking of 
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the underlying symmetry of two domain-wall states. Be¬ 
cause the loca tion of a symmetry-breaking perturbation 
is important (25), having dissipation at sites other 


than the boundary ones could still result in an exponen¬ 
tially small gap, see also recent Ref. m for a study of 
stability of edge modes to Markovian dissipation. 

We note that the scaling of the average magnetiza¬ 
tion current in the nonequilibrium steady state looks 
diffusive J 20 for this model in the gapped regime while 
higher current cumulants show anomalous non-diffusive 
scaling 153 ■ Non-diffusive scaling g ~ 1 /L 3 of the gap 
observed here (irrespective of A) is perhaps an additional 
indication speaking in favor of anomalous transport prop¬ 
erties of the gapped Heisenberg model. 
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E. Fast boundary-driven relaxation 

So-far in all models studied the gap was ~ 1/L 3 or 
smaller. On the other hand the Lieb-Robinson argument 
puts a larger upper bound on the gap of ~ 1/L. One 
can wonder whether this upper bound can be saturated? 
We are going to demonstrate that there are models with 
g ~ 1/L. To achieve that we are going to take chaotic 
models with boundary driving. 


1. Staggered XXZ with boundary dephasing 

We are first going to consider the XXZ chain in a stag¬ 
gered magnetic field, 


u = X 

3=1 


+ <T J CT J+1 


^* 5+1 


3=1 


(26) 


with staggered field having a period of 3 sites, bj = 
(— 1, — 0, — 1, — 0,...), for which the Hamiltonian is 

quantum chaotic |58j and shows diffusive magnetization 
transport [55i. Dissipation will be dephasing on the 1st 
and the last site. Similarly as in other models, z and r are 
conserved and one can look at the gap in each r-particle 
sector with z = 0 . 

In the 1-particle sector the staggered field has no in¬ 
fluence on the asymptotic gap and Eq.(22), with scaling 
~ 1/L 3 describes asymptotic g well, see Fig. 
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In the half-filling sector we are first going to evaluate 
perturbation theory for small 7 . The procedure is exactly 
the same as for the XXZ model without the field - one 
has to calculate the largest non-zero eigenvalue C\ of the 
matrix R written in Eq. ( f20j ) using eigenstates of the 
staggered XXZ chain - with the gap then being given by 
g ~ 7 C 1 . Results are in Fig. [17] . For A < 1 eigenvalue 
Ci scales as ~ 1/L, while for A > 1 (data not shown) it 
is exponentially small in L. From Fig. we also see 
that the validity of small -7 approximation shrinks with 
growing L. 

For non-small 7 we numerically calculated gaps in the 
half-filling sector, see Fig. [18]. Besides exact diagonaliza- 


FIG. 16. The gap g for the staggered XXZ chain with bound¬ 
ary dephasing (7 = 1) in the 1-particle sector. Dashed black 
lines are Eq.(22), having asymptotic decay ~ 1/L 3 . 
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FIG. 17. (Color online) Staggered XXZ chain with A = 0.5 
and weak boundary dephasing. (a) The largest non-zero 
eigenvalue of R, Eq. (20 1 , determining the gap for small 7 . 
Asymptotic behavior (full line) is tts 3.5 /L. (b) With increas¬ 
ing size the convergence radius y c upto which perturbation 
theory holds decreases. 


tion we also used tDMRG to obtain the gap, observing 
relaxation of tr(p(t)cr^ 2 cr^ 2+1 ), however, the required 
matrix dimension increases with L rapidly and one can 
not go to large system sizes. We see that for A < 1 the 
gap scales as ~ 1/L, which is different than without the 
staggered field, when it is ~ 1 /L 3 (Fig. [Til})). For A > 1 
though the gap is exponentially small, the same as with¬ 
out the staggered field. Explanation is again in terms of 
localized modes. In the gapless phase the smallest gap is 
from the 1 -particle sector, in the gapped phase it is from 
the half-filling sector. The global gap therefore scales as 
~ 1/L 3 for A < 1 and as ~ exp (— aL) in the gapped 
phase. 

The fact that the gap remains ~ 1/L 3 in the 1-particle 
sector despite chaoticity is not surprising. What is inter¬ 
esting and puzzling is that in the half-filling sector and 
for A < 1 the gap looks ~ 1/L (at least for the sizes 
available, Fig. 18), despite chaoticity and diffusive trans¬ 
port. We do not understand at present how such fast 


relaxation is compatible with diffusion, we note though 
that relaxation towards the steady state and transport 
properties of the steady state are in principle two sepa¬ 
rate properties. 
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FIG. 18. (Color online) The gap g for the staggered XXZ 
chain with boundary dephasing (7 = 1 ) in the half-filled sec¬ 
tor. The full black line in the main plot is 1 /L; in the inset 
it is oc exp (—1.5L). Full symbols are obtained by exact diag- 
onalization, empty with tDMRG. 


FIG. 20. (Color onli ne) The gap g for the chaotic XXZ model 
with staggered field (261 and boundary magnetization driving 


_p, g = 0.1, r = 1. Full symbols are exact diagonalization 

while empty are obtained from tDMRG. Straight lines suggest 
cx 1/L and ocl/L 2 asymptotic behavior. 


2. Magnetization-driven staggered XXZ 

It could be that the above fast relaxation is due to con¬ 
servation of r. We therefore take the same XXZ chain 
with staggered field as above, Eq. (26 1 , but this time 


with boundary magnetization driving of Eq. (23) instead 


of with dephasing. The value of driving is chosen to be 
g = 0.1. Now the Liouvillian conserves only z, and the 
reported gaps are for the z = 0 sector (eigenvalues in 
other sectors have larger gaps; sectors z ^ 0 also do not 
contain any steady state). The steady state is nontriv¬ 
ial and represents a nonequilibrium state with nonzero 
magnetization current. For small dissipation one can 
again use perturbation analysis in dissipation strength 
T. Everything is similar to the case of the XXZ model 
without staggered field. The perturbation matrix R is 
given in Eq. (24). For the XXZ model with a staggered 


field eigenvalues of R do dependent on g, however depen¬ 
dence of ci is very weak, with the correction at g = 0.1 
being less than 1%. Therefore one could again use a sim¬ 




FIG. 19. (Color online) Staggered XXZ chain with weak 
boundary magnetization driving, (a) The largest non-zero 
eigenvalue of R, Eq. (241, determining the gap for small T, 
scales as ~ 1/L. (b) Agreement between perturbative Ci from 
(a) with the exact gap g(r), all for A = 0.5. Convergence ra¬ 
dius F c decreases with L. 


plified expression for R given in Eq. (25). Data in Fig. 19 
show that ci, and therefore also g for small F, scales as 
~ 1/L, with the range of validity (convergence radius r c ) 
decreasing with system size L. Observe also that the ef¬ 
fect of staggered field is much more pronounced than in 
the staggered XXZ with boundary dephasing (Fig. 171. 

Going to non-small values of dissipation T, Fig.|20|, the 
asymptotic scaling of the gap seem to be ~ 1 /L for A < 
1, while for A > 1 it looks like g ~ 1/L 2 , in both cases 
though convergence is less clear than in other models. 
For A < 1 the asymptotic gap certainly seems to be 
larger than ~ 1/L 2 , showing that the Liouvillian gap for 
a boundary-only dissipation does not necessarily reflect 
diffusivity of the Hamiltonian. 


3. Tilted I sing with boundary dephasing 

The last model that we are going to study is again 
quantum chaotic one, but one that no longer conserves 
z. We take the Ising chain in a tilted magnetic field, 


L—l 

h = y. - 2 »: 

i=i 


j a j +1 


+ b x cr* + b z a^, (27) 
3 =1 


which is quantum chaotic |60j for generic field direction 
(we use b x = 3.375 and b z = 2) and with diffusive energy 
transport EM- For small dephasing strength 7 the 
gap will again be equal to g ~ 7 C 1 , where ci is the largest 
non-zero eigenvalue of matrix R (20 1 . Numerically com¬ 
puted ci for sizes L < 15 show scaling ci oc 1/L“ with 
a ss 0.85 — 1.0 (data not shown). 


We have also numerically calculated g (Fig. 21) for a 
system with boundary dephasing 0 of strength 7 = 1 . 
For smaller L we used exact full diagonalization or the 
Arpack package, while for L > 10 we used tDMRG m, 
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FIG. 21. (Color online) Gap g for the tilted Ising model (271 
with boundary dephasing of strength 7 = 1. Points are nu¬ 
merical calculations: full squares for exact diagonalization, 
and empty squares for tDMRG. 


Sector 

XX+deph. 

Bulk dissipation 

XXZ+deph. 

XX+hopp. 



A < 1 

A > 1 


1-particle: 

3(7 > 7c) 

1 /L 2 

IO 

1/L 2 

1/L° 

3(7 < 7c) 

1/L° 

1/L° 

1/L° 


7c 

1/L 

1/L 

1/L 


Half-filling: 

3(7 > 7c) 

1/L 2 

IO 

1/L 2 


3(7 < 7c) 

1 /L° 

1/L° 

S3 1/L 0 ' 8 


7c 

1/L 

1/L 12 

e-fci 



TABLE I. Asymptotic scaling of the gap g with system size L 
for studied spin chains with bulk dissipation. Gap scaling for 
perturbatively weak dissipation, 7 < 7 C , as well as for non¬ 
small dissipation is listed. Behavior of with increasing L 
is listed in two rows with heading “ 7 C ”. Approximate sign « 
means that the scaling was inferred from small L and that 
the observed scaling is perhaps not yet the asymptotic one. 


inferring g from the observed relaxation rate of to¬ 
tal energy, the initial state being a pure Neel state 
|0101.. .)(0101... |. For this model we could get g with 
tDMRG for significantly larger systems than for the other 
two chaotic models. We can see from Fig. [2l] that the 
decay seems to be algebraic 1/L“ with a being between 
1 — 2. While smaller L seem to be described by 1/L 125 , 
for larger we get nicer fit with 1/L 15 . We in principle can 
not exclude the scaling becoming ~ 1 /L 2 for still larger 
L, however, for the available data the scaling seems to be 
distinctively different than ~ 1/L 2 . 


V. SUMMARY AND CONCLUSION 

Let us briefly summarize our results and point to some 
interesting findings and open issues. We have studied 
open quantum spin chains with two types of sites at 
which dissipation acts. The first case was chains with 
bulk dissipation in which dissipation acts on all (or most) 
sites, while the second case was chains with boundary 
dissipation. Scaling of the gap for models with bulk dis¬ 
sipation is summarized in Table [i] . 

Not surprisingly, because there are no fundamental 
limitations on the gap, one also finds all different scal¬ 
ings. Quite universally, the gap scales differently for 
small dissipation strength than for non-small dissipation 
(dissipation strength is in our cases mostly dephasing 7 ). 
Systems with bulk dissipation therefore typically undergo 
a (nonequilibrium) phase transition from a phase domi¬ 
nated by H to a phase where dissipation is dominating. 
Such transition in the decay mode has been analyzed in 
detail for the XX chain with dephasing. In our models 
the critical dissipation strength 7 C at which this transi¬ 
tion happens always goes to zero in the thermodynamic 
limit. 

For boundary driven models the gap scaling is summa¬ 


rized in Table [TT] . They all comply with a general bound 
prohibiting faster than g ~ 1/L relaxation. Compared to 
bulk-dissipated cases, here the scaling seems the same for 
small and for non-small dissipation (except perhaps for 
the tilted Ising case, and the staggered XXZ model with 
magnetization driving; in both cases though finite-size 
effects could still be at play). In the 1-particle sector the 
gap is always ~ 1/L 3 due to essentially the solvability 
of the smallest nontrivial subspace (eventhough a model 
might be chaotic in larger invariant subspaces) and the 
longest eigenmodes having wavelength oc L. 

A number of transitions in the scaling of g can be 
identified. At each such transition there is a possible 
(nonequilibrium) phase transition that would be inter¬ 
esting to explore in more detail. Changing a bulk pa¬ 
rameter, like the anisotropy A or the staggered field, can 
change the scaling (e.g., from algebraic to exponential, 
or from ~ 1/L 3 to ~ 1/L). More interestingly, the scal¬ 
ing can also change by changing a boundary dissipation 
only, that is, changing only terms that have a relative 
weight 0(1/L) in the Liouvillian C. An example is the 
gapped XXZ chain for which the gap is exponentially 
small if dephasing is at the boundary, while it is alge¬ 
braic for magnetization driving at the boundary. Such 
a transition is due to a symmetry breaking of otherwise 
protected subspace. 

One finding that needs further exploration is a fast 
g ~ 1/L relaxation in chaotic models. In the present 
work we studied only the spectral gap, without detailed 
discussion of the associated eigenvector properties. Of 
particular interest is locality of decay modes and with 
it connected relaxation of local observables. Namely, it 
has been observed in Lindblad equations as well as in 
classical systems [3] that (certain) local observables can 
relax in a time that is smaller and scales differently than 
the global gap. For instance, in the XXX chain with 
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Boundary dissipation 


Sector 

XX+deph. 

XXZ+deph. 

XXZ+^i. 

stagg. XXZ+deph. 

stagg. XXZ+/I 

tilted Ising+deph. 


A < 1 A > 1 

A < 1 A > 1 

A < 1 

A > 1 

A < 1 A > 1 


1-particle: 








3(7 > 7c) 

1 /L 3 

1/L 3 1/L 3 

not conserved 

1/L 3 

1/L 3 

not conserved 

not conserved 

g( 7 < 7c) 

1 /L 3 

1/L 3 1/L 3 


1/L 3 

1/L 3 



Half-filling: 








g( 7 > 7c) 

1 /L 3 

1/L 3 e- aL 

1 /L 3 1 /L 3 

1/L 

e~ ctL 

S3 1/L S3 1/L 2 

Ri 1/L 1 5 

g( 7 < 7c) 

1/L 3 

1/L 3 e~ aL 

1/L 3 1/L 3 

1/L 

e-ocL 

1/L S3 1/L 

S3 1/L 0 ' 9 

7 c or r c 

L° 

L° 

r c ~ l° 

7c -)■ 0 

S3 L° 

r c -> 0 r c -)• 0 



TABLE II. Asymptotic scaling of the gap g with system size L for studied spin chains with boundary dissipation. For systems 
that do not conserve particle number r the global gap is listed under the “half-filling” heading. When known, we also list 
behavior with increasing system size of the convergence radius y c (or F c for “magnetization” driving) of perturbation series in 
dissipation. “Deph” in the model description denotes dephasing dissipation, Eq. 0 . “p” magnetization driving, Eq. (|23|). 


Lindblad magnetization driving the gap scales as g ~ 
1/L 3 whereas local magnetization and current relax as 

- I/L 3 / 2 |17j. 

Appendix A: Uniform-mixture steady state 

In all systems studied that conserve magnetization as 
well as the number of particles in the bra and ket, that 
is z and r (0), the Liouvillian eigenproblem has a block 
structure. In the z = 0 sector one has a steady state 
in each r—particle sector, and that steady state is an 
equal mixture of projectors to all basis states (an ergodic 
diagonal state), 

P= ^J2 l+K+l’ ( Al ) 

with | ipj) having r spins in state |0) and L — r in state 
|1). The two most interesting subspaces are the 1-particle 
and the L/2-particle (half-filling). 

1. One-particle sector 

The simplest case is L = 2, for which the steady state 
is 

p 2 = I(|10)(10| + |01)(01|). (A2) 

Taking a local operator basis composed of 

{|0)(0|, 11)(11, 10)(11, |1)(0|} one can “vectorize” the 

operator p, writing it as a vector in a Hilbert space of 
operators. Doing that on P 2 , it can be written as 

|p 2 » = -L(|10» + |01») J (A3) 

where we use the notation (•)) to denote vectors in 
the space of operators, written in the operator basis 


{10», |1», |2», |3» } ee {|0)(0|, |1)(1|, |0)(1|, |1)(0|}. The 
steady state in the 1—particle sector on L spins, p^, is 
simply 

\pl)) = 7 E (|1 ° '" 0)> +1010 ■'' 0>) + "' +10 ''' 01))) ’ 

(A4) 

which is the so-called IV-state. One can immediately see 
that for a bipartite cut after p spins the Schmidt decom¬ 
position is of rank 2 with the two eigenvalues (squares of 
Schmidt coefficients) being 1 — and j, irrespective of 
p. In the operator space such a steady state is therefore 
of finite rank 2 (it is weakly entangled) regardless of the 
system size L. 


2. Half-filling sector 

The half-filling sector is composed of states with half 
of the spins pointing up and half pointing down. Total 
magnetization is therefore zero. Starting again with a 
simple example for L = 4, we have the steady state 

Pa = i(|0011)<0011| + |0101)(0101| + |1001)(1001| + 

+ |0U0)(0U0| + |1010)(1010| + |1100)(1100|), (A5) 

or, written as a vector, 

|p 4 » = ^0 OO11 » + l° 101 » + I 1001 )) + l° 110 » + 

+ |1010))+ 11100))), (A6) 

where \pl)) is a uniform superposition of all ( L ^ 2 ) basis 
states with zero magnetization (for simplicity, we assume 
even L; for odd L and the largest sector one has to re¬ 
place L/2 with (L + l)/2). Regarding (operator) Schmidt 
decomposition, for a cut after p = 1 sites we see that the 
Schmidt rank is 2 with both eigenvalues being 1. For 
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p = 2, i.e., a cut of maximal size for p 4 , we have decom¬ 
position 

|p 4 » ~ |00»111)) + 111)) 100)) + (A7) 

+V4^(|01)) + |10)))^(|01)) + |10))), 

and therefore the state is of rank 3 with eigenvalues being 
g, g, We can see that the eigenvalue prefactors (1,1,4) 
are actually of a combinatorial nature, resulting from the 
number of combinations of distributing k ones on p sites, 
e.g., 1 = (o) ©’ 1 = ( 2 ) ( 0 )’ and 4 = ( 1 ) (i)' Generalizing 
to a bipartite cut of \pl)) after p sites, one has p + 1 


nonzero Schmidt coefficients with the eigenvalues being 
(k) ( L/ 2 -k ) / (l) 2 ) > k = 0,...,p. We can see that the 
largest Schmidt rank is for a half-cut, and is equal to 
L/ 2+1. Operator Schmidt rank of the steady state in the 
half-filling sector grows linearly with the system size L. 
We remind that a constant [52] or a linear [53] Schmidt 
rank is a sign of an exact solvability of a steady state. 
However, our results show that the exact solvability of the 
Lindblad steady state in general does not tell us anything 
about solvability of (closest) decay modes or the behavior 
of the gap. 
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